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Abstract 

We study the Sherrington-Kirkpatrick model, both above and below the De Almeida 
Thouless line, by using a modified version of the Parallel Tempering algorithm in which 
the system is allowed to move between different values of the magnetic field h. The 
behavior of the probability distribution of the overlap between two replicas at different 
values of the magnetic field h,Q and hi gives clear evidence for the presence of magnetic 
field chaos already for moderate system sizes, in contrast to the case of temperature 
chaos, which is not visible on system sizes that can currently be thermalized. 
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1 Introduction 



The Sherrington-Kirkpatrick (SK) model was introduced quite a long time ago [0J as a mean 
field model for spin glasses. Its proposed analytical solution || displays intriguing features 
such as an infinite number of pure states in the glassy phase, described by an order parameter 
which is the non trivial probability distribution of the overlap between two states, P{q)- After 
more than twenty years this solution is still the subject of works aiming at establishing it in 
full mathematical rigor || |J, whereas long standing open issues concern the study of the 
corrections to the mean field approximation below the upper critical dimension H and the 
very applicability of the mean field picture to short range realistic spin glasses |J . 

An interesting question concerns the way in which the states reorganize themselves when 
the system is subjected to a small perturbation Sp of an external parameter, in particular the 
temperature T or the magnetic field h. There is the intriguing possibility of p chaos, namely 
that states at p and states at p + Sp are as different as possible in the thermodynamic limit. 

The possible presence of temperature chaos in the SK and related models is an old subject 
of investigations p|-|IO|] that recently received a lot of attention both analytically and numer- 



ically |Tl]]-[l(J. From a very recent analytical computation |L7[] it turns out to be present but 
to be of the ninth order in perturbation theory, a very weak effect, extremely difficult to be 
numerically observed on the system sizes one is currently able to thermalize. 

The aim of this paper is to investigate the appearance of chaos with increasing system sizes 
(a question that cannot be addressed by existing analytical techniques, that are restricted to 
the asymptotic N — > oo regime), in a case where chaos is strong, namely the case of magnetic 



field chaos. The presence of magnetic field chaos was predicted already twenty years ago |L8 



(see also [0, [Hi). From the numerical point of view, it was observed in a previous work |§ 
from a study of the behavior of the second moment of the probability distribution of the 
overlap Ph ,hi(q) between replicas at h = and hi ^ 0. This pioneering paper can however 
be criticized since many data points are on the wrong side of the De Almeida Thouless (AT) 



line [p!9[ . We will revisit the problem by looking (in the SG phase) at the distribution PhoMil) 



itself, a quantity whose interpretation is simpler than the moments. 

More in detail, in terms of the probability distribution of the overlap between two replicas 
at different values of the external parameter h Q and hi = h + 8h, chaos has a very clean 
signature. Taking for simplicity the case h Q = 0, for small volumes Po,sh{q) has two peaks, 
and is very similar to Po,o{q) on the same volume. As the volume grows, a peak develops 
around the minimal value of the overlap q m = 0, in such a way that for very large volumes 
Po,o(q) ~ S(q). In the temperature chaos case, this chaotic peak is hardly visible with current 
computers and algorithms. Our aim is to determine if and how this "chaos peak" scenario 
takes place in the case of h chaos, that is believed to be much stronger than T chaos. 

To this aim, we perform numerical simulations of the SK model at T = 0.6T C , both above 
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and below the AT line, by using a modified version of the Parallel Tempering algorithm [2D. [H] 
in which the system is allowed to move between different h values at fixed temperature. 



2 Model and Observables 

The Sherrington-Kirkpatrick spin glass model ^3j is described by the Hamiltonian 

Hj = JijCTiCTj-k ^ CTj, (1) 

l<i<j<N l<i<N 

where <7j = ±1 are Ising spins, the sum runs over all pairs of spins and are quenched 
identically distributed independent random variables with mean value Jy = and variance 
l/N. We take J {j = ±N~ 1 / 2 . 

In order to measure the probability distribution of the overlap P(q) one usually considers 
two independent replicas {<Ji} and {r,} evolving contemporaneously and independently (at 
the same temperature and at the same value of the magnetic field): 

1 N 

JV i=l 



P(q) = Pj(q) = (5(q-Q)), (3) 

where the thermal average (•) corresponds to the average over Monte Carlo time in the 
simulation whereas (■) stands for the average over the Jij realizations. This is the order 
parameter in the glassy phase, which in the thermodynamic limit behaves as 

5(q-q EA ) \h\>h AT (T) 
P(q) = { x m 5(q-q m ) + P(q) + x M 5(q-q EA ) < \h\ < h AT {T) (4) 

| \P(q) + P(-q)] + x -f [5(q - q EA ) + 5(q + q EA )\ h = 0, T < T c 



where Ii A t{T) is the critical value of the magnetic field signaling the AT line, with h A x(T) ~ 
(4/3) 1/2 (T c - T) 3 / 2 for T -> T~ (T c = 1 in this model) 0. In the glassy phase, the stable 
solution corresponds to a full replica symmetry breaking (FRSB), i.e. to a non-trivial P(q) 
with a continuous distribution P(q) between two (^-functions at values q EA and q m respectively. 
For T — > T~ one finds that x m oc q m oc h 2 ^ 3 , (q EA — q m ) oc {%m ~ x m) oc (}i A t(T) — h). Note 
that at h = the function P(q) is symmetric, reflecting the symmetry of the system for 
{o"j} — > {— cTj}, and the 5-function in q m disappears. 

The interesting quantity to study when looking for chaos is the probability distribution 
of the overlap between two replicas which evolve at different values of the magnetic field, ho 
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and hi = h + 8h, definable as 



Ph ,hM = (<*(?- Qh ,hx))- 



(5) 



It is expected to become a 5-function in the thermodynamic limit, where in presence of chaos 
states are as different as possible and accordingly their mutual overlap approach the minimum 
possible value, i.e. q m (ho) (which is zero for ho — 0). This happens certainly in the iV — > oo 
limit as soon as the condition (hi — ho) 2 N » 1 is verified |18|. 

In finite dimensions, one can define the overlap correlation function Ch ,h 1 {\^i — r j|) = 
( cr i <J j)( r i r j) which decays exponentially with a correlation length that was evaluated (in d > 
8, i.e. above the upper critical dimension of the model) Q to be 6i =o,/n ^ hi~ 2 ^ 3 and 
Ch ^oM=ho+5h oc h ~ 1/6 {Sh)~ 1/2 respectively. 

Adimensional ratios of momenta, such as 



A 2n (h ,hi,T) 



B 2n (h ,h u T) 



((q-(q) 



ho,hi 



2n, 



> 



\ 2n / \ 2n 



\ 2n 



< ?-<9>i 



2n. 



(6) 



(7) 



have been introduced in order to look for chaos in reference ||, where it was argued that 
they should scale as f(Nh-y ) for h = and approach zero for iV — ► oo, namely that there 
is magnetic field chaos. 



The finite size corrections to the asymptotic behavior of PhoMil) were computed in [10 
by considering two replicas, at different values of the magnetic field, constrained to have a 
fixed overlap q. The constraint causes a free energy excess for q ^ q m given by Af — f(q — 
q m + Sq) - f(q = q m ), with 



A/ 



f /2187V/ 3 5q 2 hl /3 

5q 2 h 5h 
V2 



hn = 



ho 7^ 0, 5h = \hi — ho\ « h 



(8) 



Correspondingly one has Ph ,hi(q) ^ ex P( — ^A/), i.e. a Gaussian with variance (q 2 )oM oc 
(Nh 8 / 3 )- 1 for h = 0, in agreement with the above scaling law. 
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3 Parallel Tempering in Magnetic Field 



The Parallel Tempering (PT) or Multiple Markov Chain Method is a widely used numerical 
algorithm particularly efficient for simulating (some) systems with a corrugated free energy 
landscape. The basic idea is that the system at equilibrium, instead of being trapped in a 
single low temperature valley is allowed to move at higher temperatures where the landscape 
is trivial and to return at low T in a different valley. This can be achieved by considering n 
replicas of the system, each one at a different temperature in a given set (of temperatures), 
and by allowing exchanges of temperatures between nearest neighbor replicas with the usual 
Monte Carlo probability. 

Here we consider a set of replicas at different values of the magnetic field, both above and 
below the AT line, allowing exchange of h values between nearest neighbor replicas with the 
appropriate probability 

P ({h u {a 1 }; h 2 , {a 2 }} - {h 2 , {a 1 }; h u {a 2 }}) = j \^ Htot ^ > ° (9) 



where 



and therefore 



#« = E E Jvofof-h^ e < ( 10 ) 

o=l l<i<j<N l<i<N 



( N N \ 

E^-E-n 

u=l i=l / 



A2T« = -fa - fc) E^-E< (11) 



In principle this h-PT method should be efficient for thermalization like the usual T-PT 
method, since the landscape is trivial above the AT line. It is moreover a well suited method 
for the kind of numerical study we are interested in, since one can easily measure PhoMio) by 
considering two or more independent sets of replicas. However, as we are going to discuss in 
detail, we find that its efficiency rapidly decreases when simulating large system sizes. 

We studied the cases N = 64, 256 and 1024, taking a set of n = 49 equally spaced 
magnetic field values between h = ±\h max \ = ±0.6 at the temperature T = 0.6, where the 



AT line occurs at the critical value Hat{T = 0.6) ~ 0.382 [24 . 

We alternate one sweep of each replica with the usual Metropolis algorithm and one sweep 
with the PT algorithm. 

The probability for two replicas to exchange their magnetic fields is related to the overlap 
between the corresponding histograms of the magnetization P(m) that we check to be large 
enough (see [Fig. 1]) even for N = 1024. However some single sample Pj(m) display two 
peaks at ±m 7^ when h — (see [Fig 1]). As a result, replicas can separate into two distinct 
subsets, one evolving in the phase space with positive and the other with negative magnetic 
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Figure 1: On the left we plot the disorder averaged probability distribution of the magneti- 
zation P(m) at the 21 different central h values of the set (i.e. from h = —0.25 to h = 0.25) 
for the largest considered system sizes N = 1024. On the right we present Pj(m) at h = 
for a two-peak sample for N = 1024 again. In this last case the errors are roughly evaluated 
as the difference between the values measured in the second quarter and in the second half of 
the run. 



field values (the probability for a replica which arrive at h = with m ~ -m < to move at 
a positive 5h value is of order exp(— /3mo<5/iiV), much smaller than the usual exp(— (3x$h N)). 
This happens for some N = 1024 samples. In order to avoid such a problem, we add a new 
possible global movement, allowing a replica at h — to reverse the sign of all its spins with 
probability 1/2. 

Each run is divided into two equal parts and we check thermalization by comparing the 
data obtained in the second part with the ones of the second quarter, looking in particular at 
the behavior of P^hM)- We perform 50.000+50.000, 100.000+100.000 and 300.000+300.000 
h-PT steps for N = 64, 256 and 1024 respectively. In the N = 1024 case we also performed 
independent runs with temperature PT for 64 disorder samples at h = and h = 0.3, 
obtaining indistinguishable results for P(q). 

We simulated four sets of replicas evolving contemporaneously and independently (i.e. 
49 x 4 = 196 replicas). Data are averaged over 256 disorder configurations for each system 
size, and statistical errors are evaluated from sample-to-sample fluctuations by using the 
Jack-knife method. The program was multi-spin coded with 64 different sites of the system 
in the same computer word and the whole simulations took about 5500 CPU hours (in the 
largest part used for N = 1024), i.e. about one week when running over 32 processors on the 
COMPAQ SC270 (the program can be easily parallelized by running different samples over 
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different processors). 

In the N = 64 and 256 cases the algorithm works quite nicely, as can be seen from the 
number of tunnelings, namely the number of times that each replica moves from one extrema 
of the set (of h values) to the other and back, which is about M = 15 4- 20 (in the second 
half of the run). On the other hand, already for N = 1024, despite the 300.000 PT steps of 
the second part of the run, this number drops to M = 5 -4 6 and in nearly 1/4 of the samples 
there is at least one replica which is unable to go from h max to h m i n and back in the whole 
considered interval (in a few cases most replicas never did it). 

In the case of temperature PT, the corresponding (average) number of tunnelings are 3780, 
1590 and 455 respectively (in runs of 400.000 steps starting from equilibrium configurations, 
with a set of 38 temperatures between T max = 1.325 and T min = 0.4 at h = 0.3 for N = 64, 256, 
and 1024). Clearly the number of tunnelings decreases much faster with system size in the 
h-PT case. This is presumably linked to the early appearance of magnetic field chaos. As 
we will discuss in detail in the next section, Pho^io) starts to approach a 5-function, i.e. 
its thermodynamic limit behavior, already for magnetic field differences of order 0.15, for 
N = 1024. This means that the corresponding phase spaces are very different and that 
an algorithm based on global movements between different values of h can not work well. 
Similarly the efficiency of temperature PT should drop down for extremely large systems due 
to temperature chaos finally coming out. 

The bottom line is that N = 1024 is the largest size we are able to efficiently thermalize 
at T = 0.6 by using the h-PT algorithm, to be compared with the four time larger iV = 4096 
that can be thermalized down to T = 0.4 with the temperature PT algorithm at zero magnetic 
field. 



4 Results and Discussion 



4.1 On the finite size corrections to the P(q). 

The function P(q) is shown in [Fig. 2] for h = 0.0, 0.1, 0.2 and 0.3. At h — 0.0 it agrees nicely 
with the expected behavior, whereas it is strongly affected by finite size effects for non-zero 
magnetic field. This is in qualitative agreement with the theoretical finding [^5|] that the finite 
size corrections of P(q) are an order of magnitude larger in the q < q m region than in the 
q > qEA region, namely, 



ln(P(g)) 
N 



'^m(qm — q) 3 



q « q m 
q » qEA 



;i2) 



with A m << Xea- The behavior for q > qEA was tested (for h = 0) in |26[ and |27 . 
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Here we find that for the considered sizes P(q) has a visible tail in the q < region for h as 
large as 0.3. Moreover the peak that should correspond to the thermodynamic limit S(q — q m ) 
is not visible and the expected exp(— NX m (q m — q) 3 ) behavior is swamped by the reminiscence 
of the q = —qEA peak, still clearly visible at h — 0.1. For increasing magnetic fields the weight 
of the q = q m peak should increase (and the reminiscence of the q = —qEA peak fade away) 
but q m approaches qEA making difficult to distinguish between the two peaks. This kind of 
strong finite size effects in magnetic field were already observed in finite dimensional spin 



glasses [28, . Larger system sizes and / or lower temperatures would be needed in order to 
see the correct large volume behavior. 

On the other hand, we note that in our data qEA (defined as the location of the maximum 
of P{q)) is practically independent of the field, as predicted by Parisi theory (in the infinite 
volume limit). We obtain qEA — 0.53 for N = 1024, where a recent analytical computation 



24| gives the asymptotic value qEA — 0.505 (independently of h). 



4.2 On magnetic field chaos 

In order to find evidences for magnetic field chaos we analyze the behavior of Ph ,hi(^)- We 
first consider the case h = 0.0 (then P ,hi(<7) is still symmetric for q — > —q) and let /iitake 
the values 0.1, 0.15, 0.2 and 0.3 (see [Fig. 3]). Already for h\ = 0.15 we find clear evidences 
for a chaotic behavior when looking at the N = 1024 data. This is very different from the 
situation one finds when looking for temperature chaos [O] where Pt ,Ti(<7 ~ 0) does not 
show a clear peak corresponding to the thermodynamic limit 5(q) for (7\ — T ) as large as 0.2 
and sizes as large as N = 4096. It is remarkable that the appearance of magnetic field chaos 
with increasing system sizes is a very sudden phenomenon: chaos is elusive for N = 256 and 
blatant for N = 1024. 

On the other hand, to get a nearly Gaussian behavior we have to consider at least N = 1024 
and h\ values as large as 0.3, but the variance is more than an order of magnitude larger than 
the one predicted by (||). Our data suggest that the support of P ,/i shrinks to zero as iV 
grows, and the chaotic q ~ peak dominates more and more the distribution. 

Moreover, we find that A 2n (ho, hi) and B 2n (ho, hi), which decrease for increasing sizes as 
soon as hi > 0, are in agreement with the expected scaling law 0, i.e. f(Nh^ 3 ). We consider 
in particular 

#fr M ,T)Jk^±hs±, (13) 

<(9-<9>fa,ta) >, 



h ,h 



which is plotted in scaling form in [Fig. 4]. In the limit l/(Nh^ 3 ) << 1 the scaling function 
is approaching the asymptotic regime f(Nh^ 3 ) oc l/(Nh^ 3 ) in qualitative agreement with 
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Figure 3: The probability distribution of the overlap Pho,hi(o) between replicas evolving at 
different magnetic field values, with ho = 0.0 and h\ = 0.1,0.15,0.2 and 0.3 respectively, for 
the considered system sizes. 
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the first-order perturbative result This shows that the asymptotic regime is indeed 

approached in our data, and that we can safely deduce that limA^oo B 2 (0, hi ^ 0, T) = 0. 
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Figure 4: On top, the behavior of B2(h , hi) for ho = 0.0 as function of h\ (left) and as function 
of the scaling variable l/(Nhi ) compared with the asymptotic behavior oc l/(Nhi ) for 
l/(Nh^ 3 ) << 1 (log-log plot on the right). On the bottom, the behavior of PhoM^) f° r 
ho = 0.0 as function of hi (left) and as function of Nhi compared with the asymptotic 
behavior oc J Nh^ 3 for Nh^ 3 >> 1 (log-log plot on the right). 



We conclude the analysis of the ho — case by looking at Ph ,hi(®) as a function of hi. 
It increases when considering increasing sizes (apart from the very small h values, where 
there are clearly strong finite size effects) and scales roug hly as f(Nhl /3 ) (see [Fig. 4]) with 
f(Nh^ 3 ) approaching the expected behavior oc \J Nh^ 3 for Nh^ 3 >> 1. We also note that 
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Figure 5: The probability distribution of the overlap Ph ,hi{q) between two replicas evolving 
at different magnetic field values, with h = 0.1 and hi = 0.15,0.2,0.25 and 0.3 respectively, 
for the considered system sizes. 
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The behavior of B2(h , hi) for h = 0.1 as function of hi (left) and as function of 
- h ) 4 ) compared with the asymptotic behavior oc l/(N(hi — ho) A ) (log-log plot on 



though we have plotted data only for hi < 0.4 these scaling laws appear satisfied also when 
including data corresponding to hi values on the other side of the AT line. 

Next we consider h Q = 0.1, h x = 0.15,0.2,0.25 and 0.3 (see [Fig. 5]). The P hoM (q) 
is still expected to approach a 5-function in the thermodynamic limit, now centered in q m 
(q m (h = 0.1) = 0.21 independent of T from a recent analytical study |]2"4|). However, we 
have already noted that the peak in q m is not evident in our data for (the usual) P(q) and 
correspondingly there is no clear evidence for chaotic behavior in P^^om^) • Also for the 
largest size considered, i.e. N = 1024, though a small second peak in q rs 0.05 is appearing, 
the dominant contribution is still coming from the reminiscence of the peak in qEA, whose 
mean value and height are slowly decreasing for increasing hi. As a matter of fact, when 
going to the other side of the AT line, i.e. hi > 0.4, it is this peak that survives, becoming 
centered on a definitely lower q value ((q) ho=0 1 ^ 1=0 6 — 0-18 f° r N = 1024, smaller than q m ). 

It is clear that one should look at larger iV's to get evidences for the expected Gaus- 
sian shape oc exp(— (q — q rn ) 2 /2a 2 h ) (with l/a 2 h = \^2Nh \hi — h \) in the spin glass phase. 
Therefore, it is not surprising that a quantity such as B 2 (h , hi) does not scale as a function 
of N(h\ — ho). In the case we are considering of h = 0.1 a form B2 ~ f(N(h x — h ) a ) 
still roughly works, with a ~ 4. Nevertheless, data presented in [Fig. 6] show that even for 
iV = 1024 and (hi — ho) ~ 0.3 we are very far from an asymptotic regime f(x) oc 1/x. When 
looking at larger ho, B2 definitely does not scale, for instance already at ho = 0.2 we get 
B2(N = 1024) > B2(N = 256) in the whole interval 0.2 < h x < 0.4. 
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5 Conclusions 



We performed numerical simulations of the SK model in a magnetic field at temperature 
T = 0.6 both in the glassy phase and above the AT line. We used a modified version of the 
PT algorithm in which the system is allowed to move between a chosen set of magnetic field 
values, an algorithm well suited for our purpose. We found that N = 1024 is the largest size 
one is able to efficiently thermalize with this method and we argue that this is related to the 
appearance of magnetic field chaos at this scale. 

The function P(q) shows strong finite size corrections for h > 0, with a long tail in the 
q < region that slowly disappears for increasing sizes, whereas the peak corresponding to 
the thermodynamic limit S(q — q m ) is not yet visible. 

Our main result is on the behavior of Ph ,hi (?) > which in the case of h = 0.0 shows evidence 
for chaos already at h\ = 0.15 when considering the still relatively small size N = 1024. This is 



at variance with the situation one finds when looking for temperature chaos [11|, in agreement 



with the very recent analytical finding [17] that temperature chaos is a much weaker effect. 
The appearance of the third peak in q = is accompanied by a shrinking of the support of 
the distribution. 

The expected scaling law || is well satisfied and for large Nh^ 3 P h =o,fci (?) approaches a 
Gaussian with variance oc l/(Nh^ 3 ), in qualitative agreement with the result of a first order 
perturbative computation JTDJ. 

On the other hand, when looking at the chaotic behavior for ho ^ we found ourselves to 
be still very far from the expected asymptotic regime. This is to be related to the presence 
of strong finite size effects observable also on P(q) itself. 
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